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Abstract 



In a recent paper, Newman Q surveys the literature on power law spectra in 
evolution, self-organised criticality and presents a model of his own to arrive 
at a conclusion that self-organised criticality is not necessary for evolution. 
Not only did he miss a key model (Ecolab |^J^ ) that has a clear self-organised 
critical mechanism, but also Newman's model exhibits the same mechanism 
that gives rise to power law behaviour as does Ecolab. Newman's model is, 
in fact, a "mean field" approximation of a self-organised critical system. 

In this paper, I have also implemented Newman's model using the Eco- 
lab software, removing the restriction that the number of species remains 
constant. It turns out that the requirement of constant species number is 
non-trivial, leading to a global coupling between species that is similar in 
effect to the species interactions seen in Ecolab. In fact, the model must self- 
organise to a state where the long time average of speciations balances that 
of the extinctions, otherwise the system either collapses or explodes. 

In view of this, Newman's model does not provide the hoped-for counter 
example to the presence of self-organised criticality in evolution, but does 
provide a simple, almost analytic model that can used to understand more 
intricate models such as Ecolab. 
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1 



Typeset using REVTgX 



I. INTRODUCTION 



Over the last five years, the notion that Biological Evolution is a self-organised critical 
phenomenon has gained currency, and in particular, has been championed by Bak |Q and 
Kauffman [^. Self-organised critical phenomena are characterised by a frustration between 
two processes. The archetypical example is that of a sandpile, where the process of adding 
sand to a sand pile makes the slope of that pile steeper is opposed by the instability of 
the sandpile which works to make the sandpile flatter once the slope passes a critical angle. 
One of the the most obvious manifestations of criticality is a power law spectral behaviour, 
although criticality is by no means necessary for this power law behaviour to be manifest. 

In a recent paper, Newman surveys the field to conclude that the mechanism by which 
ecosystems are driven to criticality is not well understood, but that the evidence in the fossil 
record for power law spectra of extinction event size and species lifetimes is good. Sole 
et. al. P present the best evidence yet that these distributions are power laws. However, 
Newman missed an important model of evolution, Ecolab [0,0, that is more general than 
those surveyed, and gives us the best idea yet of how evolution could be a self-organised 
critical phenomenon. 

Newman goes further to introduce his own model of evolution to make the point that 
the coevolutionary avalanches that all the other models (including Ecolab) exhibit are not 
necessary for the observed power law behaviour. He further claims that his model is not 
critically self-organised. However, the mechanism that leads to power law behaviour in 
Newman's model is precisely the same as that in Ecolab, and that mechanism is of the 
nature of a frustration between two processes that characterises Bak's sandpile model. 

II. ECOLAB 

In this section, we consider a model of evolution called Ecolab. Ecolab (perhaps unfor- 
tunately) is both the name of a model and a simulation system written by the author to 
implement that model. The ecology is described by a generalised Lotka-Volterra equation, 
which is perhaps the simplest ecological model to use. 

■hi = TiUi + ^ PijUiUj (1) 

i=i 

Here r is the difference between the birth rate and death rate for each species, in the absence 
of competition or symbiosis. f3 is the interaction term between species, with the diagonal 
terms referring to the species' self limitation, which is related in a simple way to the carrying 
capacity Ki for that species in the environment by Ki = —rijSu. In the literature (eg Strobeck 
0, Case IQ) the interaction terms are expressed in a normalised form, aij = —Ki/riPij, and 
an = 1 by definition, n is the species density. 

These equations are simulated on a simulator called Ecolab. The vectors n and r are 
stored as dynamic arrays, the size of which (i.e. the system dimension) can change in time. 
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A. Linear Stability Analysis 



Linear analysis starts with the fixed point of equation (|I|) 

-/3-V, 



n 



(2) 



where n = 0. There is precisely one fixed point in the interior of the space of population 
densities (i.e. n such that rii > 0) provided that all components of n are positive, giving 
rise to the following inequalities: 



(/3-V)^>0, Vz 



(3) 



This interior space is denoted mathematically. 

There may also be fixed points on the boundary of IR+'^'', where one or more components 
of n are zero (corresponding to an extinct species). This is because the subecology with the 
living species only (i.e. with the extinct species removed) is equivalent to the full system. 

The stability of this point is related to the negative definiteness of derivative of n at n. 
The components of the derivative are given by 



6i 



(4) 



Substituting eq (El) gives 



dfii 



(5) 



n 



Stability of the fixed point requires that this matrix should be negative definite. Since 
the (/3"V). are all negative by virtue of (^, this is equivalent to (3 being negative definite, 
or equivalently, that its eigenvalues all have negative real part. Taken together with the 
inequalities (^, this implies that 2nsp inequalities must be satisfied for the fixed point to be 
stable. This point was made by Strobeck ||^, in a slightly different form. (Note that Strobeck 
implicitly assumes that Vihi/Ki > 0, so comes to the conclusion that 2nsp — 1 conditions 
are required.) If one were to randomly pick coefficients for a Lotka-Volterra system, then 
it has a probability of 4~^^p of being stable, i.e. one expects ecosystems to become more 



unstable as the number of species increases |10 . 



B. Permanence 



Whilst stability is a nice mathematical property, it has rather less relevance when it comes 
to real ecologies. For example the traditional predator-prey system studied by Lotka and 
Volterra has a limit cycle. The fixed point is decidedly unstable, yet the ecology is permanent 
in the sense that both species' densities are larger than some threshhold value for all time. 
Hofbauer et al. |11| and Law and Blackford |jl2|] discuss the concept of permanence in Lotka- 



Volterra systems, which is the property that there is a compact absorbing set C IR 



n, 



sp 



i.e 
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once a trajectory of the system has entered Ai, it remains in AA. They derive a sufficient 
condition for permanence due to Jansen of the form: 

^PifiifiB) = ^Piin - ^^PijnBj) > 0, 3pi > (6) 

i i j 

for every fiB equihbrium points lying on the boundary (n^j = 3i), provided the system is 
bounded (or equivalently dissipative) This condition is more general than stability of the 
equilibrium — the latter condition implies that a local neighbourhood of the equilibrium is 
an absorbing set. Also, the averaging property of Lotka-Volterra systems implies that the 
equilibrium must lie in the positive cone IR.,.''''. So @ must still hold for permanence. 

Consider the boundary points fiB that are missing a single species i. Then Jansen's 
condition for these boundary points is 

ri -^Pij^Bj > 0. (7) 
j 

This set of conditions is linearly independent. Let the number of such boundary points be 
denoted by < ngp. Then the set of conditions (P) will have rank ub < v < n^p (the 
number of linearly independent conditions, so the system has at most probability 2~^'^v~^ of 
satisfying Jansen's permanence condition if the coefficients are chosen uniformly at random. 
As stability is also sufficient for permanence, the probability lies between 4~^**p and 2~'^^p~'^. 

Another rather important property is resistance to invasion. [§] Consider a boundary 
equilibrium n^. If it is proof against invasion from the missing species, then the full system 
cannot be permanent. For the boundary points that miss a single species, this implies that 
condition (0) is necessarily satisfied for permanence, along with (^. The probability of 
permanence is then bounded above by 2"^*^?"""^. 

The important point to take away from this section is that whilst a randomly selected 
ecology is more likely to be permanent than to have a stable equilibrium, the likelihood 
decreases exponentially with increase in species number. 

C. Mutation 

Adding mutation involves adding an additional operator to equation (|l]) 

h = r*n + n* /3n + mutate(/i, r, n) (8) 

where * refers to elementwise multiplication. This operator extends the dimension of the 
whole system, so is rather unusual. It is not germane to the present argument what the 
precise form of mutate is, the interested reader is referred to the previous publications 
describing it Suffice it to say, that it adds new species according to a stochastic 

mechanism, and that we would expect that the criticality result to be robust with respect 
to changes of mutation algorithm employed. 



Boundedness is ensured in this model by choosing the Pij such that Pij + (3ji < 0, Vi, j. This 
precludes symbiosis, but does allow for unstable behaviour. See for a discussion of boundedness 
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D. Self Organised Criticality 



Lets consider what happens to the largest eigenvalue of (3. Suppose initially, the system 
has a stable equilibrium, in which case all the eigenvalues have negative real part. As 
mutations are added to the system, the largest eigenvalue will increase towards zero. As 
it passes zero, the system destabilises, and the system will start to exhibit limit cycles or 
chaotic behaviour. As further mutations are added to the system, permanence is no longer 
satisfied, and an extinction event will occur. This will restore permanency to the system, 
and possibly even stability. So we have two frustrated processes opposed to each other, the 
first, mutation, which builds up ecosystem complexity, and the second being the trend to 
impermanency as ecosystem become more complex. This is analogous to the sand being 
added to the top of the pile, and the stability of the sandpile slope in Bak's sandpile model. 



III. THE NEWMAN MODEL 



Newman has presented his model of evolution in a number of papers |14,TqJlI], and is 



largely equivalent to an earthquake model presented in |T6|,|T^. In the biological context, 
the model has a fixed number of species, all of which feel an environmental stress, denoted 
by fjit), which is random variate with distribution Pstrcss(^)- Each species has an individual 
threshold Xj, such that if r] > Xi, species i becomes extinct. These extinct species are 
then replaced by new species, with thresholds randomly assigned from some distribution 
Pthrcsh(a^)- There is one further twist to the model, in that the threshold values are allowed 
to drift over time in order to prevent the model stagnating with every species having the 
maximum threshold. 

The Ecolab software allows us to build a variant of this model that allows the number 
of species to vary over time. When the model was first implemented, the system underwent 
a "mutation catastrophe" , in which the number of species exploded. This is similar to what 
happens in the Ecolab model when the mutation rate is set too high. Normally, one would 
expect that the number of speciation events should be proportional to the number of species. 
However, this leads to an excess of speciation over extinctions. 

The resolution of this conundrum is to require that the stress values rj be proportional 
to the number of species, i.e. r] = nspi]', where rj' is drawn from some distribution Pstress(^')- 
The justification for making this assumption can be seen by considering a simplified model 
of Ecolab (called Ecolab — ), described in the next section. Of course, in Newman's original 
model, constant, and so his model is consistent with this modification. 



Wilke and Martinetz examined a similar model, in which they label the mutation 



rate g, and consider finite /, rather than / = as I do here. They too note the conundrum 
of exponential growth in species number, and resolve it by introducing an arbitrary logistic 
constraint. My argument is that the reason for this logistic constraint is that species must 
interact with each other, and the greater the number and strengths of these interactions, 
the greater the stresses are that are felt by the ecosystem. 

It could be argued that the raison d'etre of the Newman model is to study the effect of 
coherent extinction through exogenous causes. However, these will always give rise to stress 
distributions that are independent of species number. However, the stress distribution will 
ultimately be dominated by the term that does depend on the species number. 
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Once the stress values depend on species number, the system self-organises so that spe- 
ciations and extinctions balance on average. A trace of Ugp can be seen in Figure |l], and the 
distribution of lifetimes is seen in Figure |^. The peak in the curve at r = 10 is an artifact of 
the simulation, and should be ignored. The distribution actually has two regions, the inner 
one 10 ^ r ^ 10'^ having a power law with exponent ~ — 1, and the outer region r ^ 10^ 
having exponent ~ —2. By running the experiment at different mutation rates, the lifetime 
A at which the distribution changed from to r^^ was found to be inversely proportional 
to the mutation rate. 

In comparing the result of my variation with the original Newman model, it should be 
noted that the power law exponent in Newman's original model is —1 out to a time 1//, and 
decays exponentially after that. In my version, the same power law exponent was observed 
out to 1/g, and then appears to change to a faster power law decay, although the error 
bars are sufficiently large not to rule out an exponential decay. In each of these models, the 
lifetime 1/f or 1/g respectively is roughly the lifetime that a maximally fit organism (one 
with a maximal value Xi) can survive before sucumbing to mutation pressures. 

IV. THE ECOLAB— MODEL 

In this section, we will consider a simplification of the Ecolab model where the interaction 
terms /^ij'^j cire replaced by a random variate r]i{t) from a suitable distribution: 

rii = {ri - rii)ni (9) 

Since rji is effectively a sum of a large number of independent quantities, its distribution will 
tend to be normal, and the deviation (controlling how large r]i gets) will be proportional to 
nsp, the connectance (proportion of nonzero elements in (S) and the interaction strength. This 
is why stresses in the Newman model must be proportional to rigp. When r]i exceeds for 
any significant period of time, species i becomes extinct. Since rii{t) is a continuous function 
of n{t) which is itself a continuous function of t, there will be a correlation T]{t)ri{t + r) > 
0, Vr < To, 3ro > 0. Equation (^) connects the full Ecolab model with the Newman model. 

In order to make the analysis simpler, we assume that Ui are real valued, rather than 
integers as in Ecolab. In order to detect when extinction happens, we take an arbitrary 
threshold a, such that if rii < a, species i is extinct. 

V. DISTRIBUTION OF SPECIES LIFETIMES 

Figure ^ shows the distribution of species lifetimes (time from speciation to extinction) 
in the augmented Newman model. This figure is not normalised, as a power law has 
an infinite integral. So the abcissa of the graph is not significant, but the slope is. The 
lines are fitted by linear regression. Authors often quote a correlation coefficient, however 
this is generally meaningless on a log- log plot. Even the value of the slope is meant to be 
an indication only, as the large relative error at high lifetime values can lead to significant 
errors in the computed slope. 

Figure |^ shows the lifetime distribution for Ecolab which has a slope of —2 for lifetimes 
less than 100, but —1 for larger lifetimes. At still larger times (r ^ 0.1/yu), the distribution 
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turns over, decaying exponentially. Previously published versions of this graph only show 
the smaller lifetime behaviour. 

Consider now the probability p{> t\x) that a species with threshold x will become extinct 
after time t = r in the Newman model. Since time is discrete in this model, this is simply 
the probability that the stress rj does not exceed x for the first r steps: 



p(> t\x) 



Pstrcss(^?)c?'7 



(10) 



Now the distribution p{> r) of species having lifetimes r is just the above quantity, 
integrated over the distribution of thresholds: 



p{> t) = p[x)p{> T\x)dx 



Pstress{v)dr] 



, . ^dx 

/ Pthreshlx)^ —d^ 
Jo "(^ 



(11) 



where ^ = Jq Pstvess{v)d'n 

Assume the following inequalities hold: 



Pthresh(a;) < KiPstressix), Vx 

> KoPstress{x) VX < X^, 3Xc 



(12) 



Without loss of generality, pthresh(a;) is taken to be the uniform distribution between and 
1, and is zero outside this interval. Pstressix) is positive for all positive x, with the large x 
tail needed to establish power law behaviour ||15|. In this case, the constants Kq and Ki 
correspond to the inverses of the maximum and minimum of pthresh(a^) over the unit interval, 
and Xc = 1. Let us introduce = Jq" Pstress{x)dx as being the change of variable equivalent 
oi Xc- In the case of uniform threshold distribution, and monotonic stress distribution, 1 — 
is the proportion of stress events that overwhelm the hardiest of species. The inverse of this 
proportion is a time scale, above which the lifetime distribution must decay exponentially. 
In order to observe power law behaviour, the stress distribution must be chosen so that 



Substituting eq ([1^ ) into (PD generates the following inequality: 



dx dx 

Ko / Pstress{x)C -Tjd^ < P{> t) < Ki / Pstress{x)C —d^ 

Jo "s Jo 

Ko^<p{> r)<K,- 



di 



r+ 1 



T + r 



(13) 



since Pstress(a;) = d^/dx and where = j^" Pstvess{,x)dx. 

Now p{t) = p(> r — 1) — p(> r), so the following inequality is obtained: 



r r + 1 



r r + 1 



(14) 
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Assuming that r ^ (1 — ^c) ^, = (1 + ^(1 — ^c) + ■ ■ ■ ~ 1, this inequahty may be simphfied: 

^ P(^) ^ ^/ , IN (15) 



r r 



1) -^v ^ - r(r + 1) 



This resuh indicates that there are two domains, the first being when r < j^^ztj^^ where 
the hfetimes distribution is a power law with exponent —2. This domain is more pronounced 
the closer Ki is to K2, ie the closer pthresh(a;) is to Pstress- The other domain occurs when 
T > j^^tk^, where any power law will have an exponent less than —1. In between, there will 
be a transition between the two domains. This result is not terribly strong, as the inequality 
can also be satisfied by any distribution falling off faster than a power law. However, it does 
contradict the results of the Time Average Approximation theory of Sneppen and Newman 



17| in the case of the Lorentzian distribution, where a power law with exponent (i.e. a 
flat distribution) is predicted. Whilst a flat distribution is manifestly rediculous, others are 
not. The TAA predicts a power law of 1/3 for a power law stress distribution with exponent 
-3/2. Figure ^ shows the observed lifetime distribution in this case, and the distribution 
never flattens out more than r~^. 

Now lets us turn our attention to the Ecolab — model to see if similar relationship can 
be derived. In what follows, the species index i is dropped. Integrating equation (^ gives 
us: 



n( 



and taking logarithms gives: 

\nn{t) = / r — r](s)ds, 







since no = 1 for all new species. 

For the species to become extinct after time t = t, we require: 



/ r — ri{s)ds > In a, \ft < t 
Jo 



(16) 



Since time is discrete in this model, 77(5) is a piecewise constant function, therefore the 
integral can be replaced by a sum so that 



t-i 



^7]i <rt-\na, \/t < t (17) 



i=0 



Now inequality ( p^Tf ) deflnes a set C W , and the probability of a species having lifetime 
greater than r if its reproduction rate is r is given by: 



p{> T\r) = Y\pstrcss{ili)driodr]i- ■ -drjr-i (18) 
Lets us flrst deal with sufficient conditions for inequality (|T7|) to be satisfied, which are: 
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< r — In (t/t, Vi < r 

< r, ascr < 1 



(19) 
(20) 



Therefore a lower bound for p(> r|r) is 

p{> r|r) > 



(21) 



Now consider the following relation: 

n(t + 1) = (1 + r - r]t)n{t) 
For the species not to go extinct before t = t, we require r/t < 1 + r, Vt < r. Therefore, 

"(r+l) 



p{> r|r) < 



Now find constants Kq and Ki so that 



(22) 



< i^iPstress(r- + 1) 



(23) 



where Pr{r) is the probability distribution of reproduction rates. Since p(> t) = f pr{r)p{> 
T\r)dr, we find: 



dr < p{> r) < Ki / Pstressl?^ + 1 



(r+1) 



^<pi>r)<^ 

T T 



dr 
(24) 



Now since p(r) = p(> r) — p{> r + 1), 



T r + 1 
{Kopl - Ki)t + Kopl 



Ki Kopl 



r(r + l) 



< P{r) < 



T r+1 
r(r + 1) 



(25) 
(26) 



Again, like the Newman model, we have two domains of power law possible, an inner 
domain where the power law is -2, and an outer domain where any power law is capped by 
-1. This is what is seen in Figure |[ 



VI. CONCLUSION 

The Newman model owes its power law behaviour to much the same mechanism as does 
Ecolab, although the assumption of constant species number hides essential interspecies 
connections. Both models demonstrate a power law exponent near —2 at small time scales. 



agreeing with the fossil record (after Sneppen et. al |]T9[); turning over into a gentler power 
law with exponent less than -1 at larger times. 
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FIG. 1. Tlgp clS Si function of time in the genralised Newman model. 
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FIG. 2. Distribution of species lifetimes in generalised Newman model with Gaussian stress 
distribution. 
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FIG. 3. Distribution of species lifetimes in Ecolab. 
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FIG. 4. Distribution of species lifetimes in 
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